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Abstract 

Using mean-field theory and high resolution Monte Carlo simulation tech- 
nique based on multi-histogram method, we have investigated the critical 
properties of an antiferromagnetic XY model on the 2D Kagome lattice, 
with single ion easy-axes anisotropy. The mean-field theory predicts second- 
order phase transition from disordered to all-in all-out state for any value of 
anisotropy for this model. However, Monte Carlo simulations result in first 
order transition for small values of anisotropy which turns to second order 
with increasing strength of anisotropy, indicating the existence of a tricritical 
point for this model. The critical exponents, obtained by finite-size scaling 
methods, show that the transition is in Ising universality class for large values 
of anisotropy, while the critical behaviour of the system deviates from 2D-</> 6 
model near the tricritical point. This suggests the possibility for existence of 
a new tricritical universality in two-dimensions. 
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I. INTRODUCTION 



The phenomenon of geometric frustration has attracted the interest of physicists due to 
the presence of degeneracy in the classical ground states arising from the arrangement of 
spins on triangular clusters [1-4]. A frustrated magnet is one in which not all interaction 
energies can be simultaneously optimized, for which the anti-ferromagnetic Ising model on 
a two-dimensional triangular lattice, is an example. The highly frustrated magnets, on the 
other hand, are the class of frustrated magnets that have an infinite number of classical 
ground states, even after removing the global symmetries of Hamiltonian. 

The classical XY anti-ferromagnet on the two-dimensional Kagome lattice constructed 
from corner-sharing triangular units and the classical Heisenberg antiferromagnet on the 3D 
pyrochlore lattice consisting of corner- sharing tetrahedra are two prototypes of the highly 
frustrated class. The discoveries, such as heavy- fermion behaviour [5], spin-ice ordering 
[6-8], spin nematics [9], spin liquid behaviours [10-12] and even novel superconductivity 
[13] in materials with magnetic sublattices of corner- sharing tetrahedra (such as spinel and 
pyrochlores), have made these structures in the focus of physicists's attention over the recent 
years. 

It has been widely accepted that no order-by-disorder mechanism can establish a long- 
range order in the Heisenberg pyrochlore anti-ferromagnet, consequently such a system re- 
mains disordered at all temperatures [14-16]. However, experimental observations have rep- 
resented an all-in all-out long-range order (consisting of four sublattices oriented along four 
[111] spin directions), for the low-temperature phase of FeF 3 in pyrochlore form [17,18]. In 
this compound, the Fe +3 ions located on a pyrochlore lattice, interact anti-ferromagnetically 
with their nearest neighbors. Since, the magnetic Fe +3 ions are in d 5 electronic configura- 
tion with a totally symmetric ground state and no net angular momentum, this system can 
be considered as a Heisenberg anti-ferromagnet and so the origin of the long-range ordered 
phase in it, has remained as a puzzle. Reimers et al have shown that, taking into account 
the interaction with farther neighbors, would cause a second order transition in this a sys- 
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tern [19]. However, they found that because of the thermal fluctuations, a co-linear spin 
ordering would be preferred rather than the all-in all-out state. Therefore, it seems that 
to stabilize a long range all-in all-out spin configuration, one should inevitably introduce 
a single-ion an-isotropic crystal field term in the model Hamiltonian. Another interesting 
aspect of the transition in pyr-Fe +3 is in its universality class. The order parameter critical 
exponent (3 has been fixed to the value 0.18(2), in neutron- diffraction experiments, which 
is nearest to the tetra-critical value f3 — 1/6 [20]. On the other hand, recent Monte Carlo 
simulations, carried on Heisenberg pyrochlore antiferromagnet with single ion anisotropy, 
have revealed the existence of a tricritical point for this system [21,22]. 

The above interesting problem motivated us to study the critical properties of its two- 
dimensional equivalent, the XY Kagome anti-ferro magnet model with single-ion anisotropy. 
The classical antiferromagnetic 0(n) models on the Kagome lattice have been studied by 
Huse and Rutenberg [23]. There, it has been shown that the Ising model (n = 1) is disor- 
dered at all temperatures, while the XY model (n = 2) represents quasi long-range order in 
a three-fold ordered parameter at zero temperature. Because the system is two-dimensional 
this quasi long-range order does not survive at finite temperatures and so transforms to 
disordered phase through a Kosterlitz-Thouless transition. The ground state of the XY 
model has the same properties as the three-state Potts model which can be mapped exactly 
onto solid-on solid (SOS) model at the roughening transition. On the other hand, the study 
of two-dimensional antiferromagnet Heisenberg model on Kagome, have been carried out by 
Ritchey et al, which resulted in a coplanar spin configurations in which there are nematic 
spin correlations with planar threefold symmetry and non-Abelian homotopy [24]. They 
have also shown that very small amounts of bond XY anisotropy are sufficient to convert a 
crossover to a topological phase transition, in which the binding of non-Abelian disclinations 
would result in a glassy behavior in the absence of extrinsic disorder. 

The Hamiltonian of nearest-neighbor XY antiferromagnet model on the Kagome lattice 
is given by: 
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if = -j5>.s,-, (i) 

(ij) 

in which J(0 and Sj denotes the unit planar vectors and (ij) indicates the nearest-neighbors. 
The ground state of this model is known to have a huge accidental degeneracy not related 
to the global symmetries of the Hamiltonian [23,25]. In any ground state of the Kagome 
lattice the spins S; acquires only three directions whose angles with respect to an arbitrary 
axis, say x-axis, differ from each other by 2n/3. Therefore, the ground state in addition 
to the continuous U(l) symmetry (due to the arbitrary simultaneous rotation of all spins) 
is characterized by a well developed discrete degeneracy of the same type as in the 3-state 
antiferromagnetic Potts model. 

The extensive degeneracy of the ground state in this model makes it extremely unstable 
towards the imposing of perturbations [26]. For instance, if one adds a single-ion easy-axis 
anisotropic term to Hamiltonian (1), all spins prefer to align along the anisotropy directions 
yielding a long-range all-in all-out state for the system. 

The goal of this paper is to determine the critical properties of an XY model on the two 
dimensional Kagome lattice with single ion easy-axes anisotropic term. For this purpose we 
employ mean-field theory and Monte Carlo simulation. 

The structure of paper is as follows. In Sec. II, we introduce a mean-field formalism 
to derive the qualitative picture of transitions in the model. Section III is dedicated to the 
Monte Carlo method based on multiple histograms and also some methods for analyzing the 
Monte Carlo data to determine the order of transitions, critical temperatures and critical 
exponents. The simulation results and discussion are given in Sec. IV and conclusion appears 
in Sec. V. 

II. MEAN-FIELD FORMALISM 

The Hamiltonian, describing the XY spins with nearest-neighbor anti-ferromagnetic in- 
teraction on a Kagome lattice subjected to single site easy-axes anisotropy, is given by : 
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H = -\ E E s? • sj - d £ E(s- • i a )\ (2) 

i,j a,b i a 

in which J < 0, D > and i,j = 1, • • •, N and a, b = 1, 2, 3 denote the Bravais lattice and 
sublattice indices, respectively. z a 's represent the unit vectors of three easy- axes directions 
in 2d plane, which are along the line connecting the corner and the center of corner-sharing 
triangular units, given by: 

7? = (0, 1) (3) 

in global Cartesian coordinates. 

To apply mean-field theory on this model, we follow the method introduced by Harris, 
Mouritson and Berlinsky [19,27]. Defining the average magnetization as Mf = (Sf) and the 
deviation from the mean magnetization as 5S® = Sf — Mf, to order 0(5S 2 ), we can write 
the Hamiltonian (Eq.(2)) as the following linear form: 

^ = ^EE M , a - M 5 + ^EE( M "-^) 2 -^EE M ?-s^-2DE(s- -nm^ a )- 

i,j a,b i a i,a j,b i,a 



(4) 



Therefore, the mean-field partition function can be written as: 
where 



(5) 



B * = J E E Mj + 2D{Mt ■ z a )z a , (6) 

jjki b^a 

in which, the summation is over the nearest neighbors. The integral in Eq.(5) can be 
evaluated easily as follows: 

J e^USI = 2tt £ e^ cos ^d9 = 27r/ (/3S»), (7) 



where Bf = |B°|. Then, assuming Kb = 1, we reach the following expression for the free 
energy: 

F = -T In Z = J -Y, E M i ■ M 5 + D E E( M " • ^) 2 - T E ln (2tt/ (^)) • (8) 

i,j a,b i a i,a 

From the mean-field free energy, obtained above, one can calculate the magnetization 
and entropy as: 

' M * = - v - f = -H^ a = - £ (io) 



For small values of B, one can expand Eq.(lO) as: 

'Bf Bf Bf 11 Bf 



2T 16T 3 96T 5 6144 T 7 
from which, by reversing the series one gets: 



(11) 



Bf = 2TMf - T(Mf ) 3 + -(Mf ) 5 + 0(M 8 ). (12) 

9 

Substituting Eq.(12) into Eq.(9) and expanding the entropy in powers of M°, enables us 
to expand the free energy as: 

F=(H}- TS 

= -ANT 1h(4tt) - - E E M ? ' M ? ~ D E( M ? ' 

i,j a,b i,a 

+ TJ2 (Vf) 2 + \{M? ) 4 - ^(Mf) 6 + 0(M 7 )) , (13) 

where we have used Eq.(4). We can also expand the free energy in terms of Fourier compo- 
nents defined by: 

M« = EM^exp(*q-Rn (14) 

q 

Jf = EE Jex P (<q • (R" - R?)) > ( 15 ) 
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where the summation in Eq.(15) is over the nearest neighbors of a selected spins. Then we 
reach the following form for the free energy per particle in terms of Fourier components: 

f(T,J,D) = F{T ' N J,D = -4Tln(4vr) 

+ \ E E M«M 6 _ q (2T5 ab - J« b ) - D ]T ]T(M q • r ) (M a q • 5 a ) 



q a 



+ 7 T EE( M ql' M q2)(M«3-M q4 ) 

a {q} 

- 4 T EE( M qi ' M^)(M^ 3 ■ M£ 4 )(M£ 5 • M",) + 0(M 7 ), 

dD «» {Q} 



where 



E = E^(Eqi)- 

{q} {q} < 

The free energy (Eq.16) can be rewritten in terms of Cartesian components of 
M q = (m^mq' 2 ) as: 

/(T, J, D) = -AT ln(47r) + ^ £ £ £(2T5 ab ^ - " D^S^m^m™ 



g ab a(3 



a a/3 {q} 



- ^ T E E E«i Q <2 )K^<4 )W ) + °( M? )> 

a a/37 {q} 



(16) 



(17) 



in which a,/3, 7 take the values 1,2. It can be seen from the above equation, that only the 
an-isotropic term D couples the different Cartesian components of M. The 2x2 matrices 
D a are given by: 
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(18) 



Thus we are left with the following coupling 6x6 matrix for the quadratic terms: 

/ & J? J? ^ 



</q= D 



J q 2 -D 2 J 23 

713 723 n3 
y J q J q 



(19) 



in which the off-diagonal matrices J q J are proportional to the 2x2 unit matrix as follows 



J q 2 = 2Jcos(|)/ 2x2 (20) 
J q 3 = 2Jcos(^|±^)/ 2x2 (21) 
4 3 = 2Jcos( ^ 2 ~ 9 - )/ 2x2 . (22) 

In deriving the above expressions, we have used Eq.(15) together with the positions of 
Kagome atoms given by their xy components. For convenience we reduce the number of 
indices (a = 1, 2, 3 and a = 1, 2) by defining a new set of indices s = 1, • • •, 6, which leads to 
a 6-component magnetization vector as: 

M q = (m^m^mJV • -mf ) = (m q ,m q ,- ■ -m q ), (23) 

from which the quadratic term in free energy can be written as: 

/( 2 ) = ^M q .4Mj. (24) 
q 

Diagonalizing the quadratic term, requires transforming to the normal modes $ q : 

< = lt (25) 
i=i 

for s = 1,2,- • -6. 

Uq is the unitary matrix that diagonalizes the coupling matrix J q , with eigenvalues A q : 

E<^? = W- ( 26 ) 

6 

in which, the unitarity condition requires: 

£t^t/^ q = <jtf. (27) 

a 

Equation (26) enables us to write the free energy as a power series in terms of normal 
modes, such that to O(0 7 ) we obtain the following expansion for the free energy: 
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f(T, J, D) = -AT ln(47r) + - £ ]T(2T - A q )0 q 0*_ q 

z g i=l 

+ jEEE ^i^^3^<i< 2 q3 <4 

s=l ijfci {q} 

+ 4 T E E E^l^^S^J^Wl^^^q5^6- (28) 
s=l ijklmn {q} 

It is clear that phase transition occurs when the sign of quadratic term of free energy 
changes. Therefore, from the above expression one finds that the spontaneously breaking 
symmetry occurs at a temperature: 

T c = imax q ,{A q }, (29) 

where max {} means the global maximum over all % and q. In the case of D = one can 
exactly diagonalize the matrix J q (Eq.(19)) and find the following eigenvalues: 

A q = -2J i = l,2 

A q = 2 J(l -y/3 + Q) i = 3,4 

A q = 2 J(l + ^/3 + g) z = 5,6, (30) 

where Q is given by: 

Q = {cos(2q x ) + cos(v / 3g K + q y ) + cos((2 - y/3)q x - q y )}, (31) 

which coincides with the result derived in Ref. [19]. The above results show that for J < 
the largest eigenvalues are degenerate and dispersionless (q-independent), such that when 
T < — J, the order parameters corresponding to all of these modes turn to be nonzero and 
we were left with a huge number of states with broken symmetry. Therefore, because of the 
extensive degeneracy of symmetry broken states, one concludes that in mean-field theory, 
no long range order can be established as the temperature decreases down to zero. The 
g-dependence of eigenvalues for D = along [1 0] direction is depicted in Fig.(l). 

For an-isotropic case (D)0) the eigenvalues of matrix J q can be obtained numerically. 
The dispersion curves for D = 0.2 and D = 1.0 along [1 0] direction has been shown in 
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Figs. (2) and (3), respectively. As can be seen from these graphs, all the degeneracies have 
been removed, so we were left with 6 distinct modes, where the highest mode has a maximum 
at q = with the value Aj = —2 J + 2D. It can be easily shown, by deriving the eigenvector 
of this mode, that this mode corresponds to all-in all-out spin configuration represented 
in Fig. (4). As a result, the mean-field theory predicts a continues phase transition from 
disordered to a long-range ordered all-in all-out state at the critical temperature T c = 
—J + D. Another interesting point is that the branch A q is independent of magnitude 
of anisotropy (D), which means that the modes describing by it, are corresponding to spin 
fluctuations perpendicular to easy- axes directions (z a , a = 1,2,3) in Hamiltonian, given by 
Eq.(2). 



III. MONTE CARLO SIMULATION 

For large values of D, spins tend to remain mainly along easy- axes directions such that the 
effective degrees of freedom flip along these axes. Therefore, one expects that the transition 
to all-in all-out state to be in 2D Ising universality class. However, when D is small, the 
transverse fluctuations normal to local easy-axes directions become larger and so this leads 
to lowering of the transition temperature as well as deviation from Ising behaviour. In this 
section we use Monte Carlo simulation, to study the phase transition of the model described 
in previous section and find the order of transitions for different values of anisotropic term 
D. 

To obtain a qualitative picture of the transitions and also the approximate location of the 
critical points, we first set some low resolution simulations. The simulations were carried out 
using standard Metropolis single spin-rotating algorithm with lattice size N = 3 x 20 x 20. 
During each simulation step, the angles of planar spins with the horizontal axes were treated 
as unconstrained, continuous variables. The random-angles rotations were adjusted in such 
a way that roughly 50% of the attempted angle rotations were accepted. To ensure thermal 
equilibrium, 100 000 Monte Carlo steps (MCSs) per spin were used for each temperature 
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and 200 000 MCS were used for data collection. The basic thermodynamic quantities of 
interest are the specific heat c = ((E 2 ) — (E) 2 ) / (NT 2 ) , the order parameter defined as 
M = | Ei, a S- • z a \/N and the susceptibility x = ((M 2 ) - (M) 2 )/(NT). 

In Figs. (5-8), temperature dependence of the energy per spin, ,the order param- 
eter, specific-heat and susceptibility have respectively been represented for J = —1.0, 
D = 0.2,0.1. As can be observed from Figures. (7) and (8), the transition for D = 0.2 
seems to be continuous, while for D = 0.1, because of sudden peaks in specific heat and 
susceptibility, it seems to be first order. However, The determination of the order of transi- 
tion requires more accurate methods, for which we will use Binder's fourth energy cumulant 
method. Once the probability density of energy (P(E,T)) is obtained, for measuring the 
thermodynamic quantities other than the energy, one can choose to work with this energy 
probability distribution and microcanonical averages of the quantities of interest. This leads 
to optimized use of computer memory. The microcanonical average of a given quantity A, 
which is a function of energy, can be calculated directly as: 

A(E) = (32) 

l^t °E t ,E 

from which, the canonical average of A can be obtained as a function of T: 

{ } ~ J2 E P(E,T) ■ (33) 

In our simulation, we use Kagome lattices with linear sizes L = 20, 24, 28, 32, 36, 40 (the 
number of sites is given by iV = 3 x L x L), such that the maximum number of spins is 4800, 
large enough for reducing the finite size effects. For each system size, at least five overlapping 
energy histograms are obtained near the transition point so that the statistical uncertainty 
in the wing of the histograms, may be suppressed by using the optimized multiple-histogram 
method [28]. This enables us to measure the location and magnitude of the extrema of the 
thermodynamic quantities with high accuracy. For each histogram we performed 5 x 10 5 
Monte Carlo steps per spin for equilibration and also 5 x 10 5 MCSs for gathering data. 
To reduce the correlation, 10 to 20 Monte Carlo sweeps were discarded between successive 
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measurements. In all simulation we fix J = — 1 and vary the value of D from 0.1 to 1.0. 
First of all, we deal with the order of transitions. 



A. Order of the transition 

To determine the order of transitions, we used Binder's fourth energy cumulant defined 

as: 

U ^ l -TTE^- < 34 > 

It has been shown that this quantity reaches a minimum at the effective transition temper- 
ature T C (L) whose size dependence is given by [29-31]: 

Umin(L) = U* + BL- d + 0(L- 2d ), (35) 

where 

U* = ^ - ( ei /e 2 - e 2 / ei ) 2 /12. (36) 

The quantities t\ and e 2 are the values of energy per site at the transition point of a first 
order phase transition and d is the spatial dimension of the system (d = 2 in our simulation). 
Hence, for the continuous transitions for which there is no latent heat (ei = e 2 ), in the limit 
of infinite system sizes, U m i n (L) tends to the value U* equal to 2/3. For the first-order 
transitions, however e\ ^ e 2 and then U* reaches a value less than 2/3 in the the limit 
L — > oo. 

The size dependences of U(L) for D = 1.0,0.18,0.15,0.13,0.1 have been exhibited in 
Fig. (9). The straight lines fitted to the data have been obtained from Eq.(35). The values 
of U* and latent heat per spin are also listed in Table. (I), from which one can see that, 
within the errors of simulation, transitions are second order for D > 0.17 and clearly first 
order for D < 0.15. The precise determination of the tricritical point is extremely difficult, 
however our results suggest the existence of a tricritical point between D/\J\ = 0.15 and 
D/\J\ = 0.17. 

12 



In the Figs. (10) and (11) the energy histograms of D — 0.2 and D = 0.1 for the size 
N = 3 x 40 x 40 have been shown, respectively. As can be seen from these figures, the 
energy histogram for D = 0.2 has one broad peak at the transition, while for D = 0.1, 
it has two well separated peaks around the transition temperature. This is in agreement 
with the results of Binder's method. Note that the small peak at the middle in Fig. (11), is 
artifact of the finite time of simulation and will vanish at large enough times. The reason 
is that at a strong first order transition point, free energy possesses two equivalent minima 
corresponding to two stable coexisting phases. For large system sizes these two minima 
are separated by a large energy barrier, so the system remains mainly around its minima 
during the time evolution, caused by thermal fluctuations, in simulation. Therefore, the 
configurations corresponding to the unstable region at the middle are rare, consequently the 
relative error for these data is large. 

As the next step we proceed to calculate the critical temperatures and critical exponents 
for continuous phase transitions, using finite-size scaling theory. 



B. Determination of T c and static critical exponents 

According to the finite-size scaling theory [32], the scaling form for various thermody- 
namic quantities such as magnetization density, susceptibility and specific heat in zero field 
are given by: 

L p/u M{tL l/v ) (37) 
X ~ IS /v K{tL 1/v ) (38) 
Coo {t) + L a ' u C{tL 1 l y ), (39) 

where t — (T — T c )/T c is the reduced temperature for a sufficiently large system at a 
temperature T close enough to the infinite lattice critical point T c , L is the linear size of the 
system and a, (3,^,6 are static critical exponents. Equations (37-39) are used to estimate 
the critical exponents. However, before dealing with the critical exponents we should first 
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determine the critical temperature accurately. 

The logarithmic derivatives of total magnetization {mL d ) are important thermodynamic 
quantities for studying critical phenomena and very useful to high accurate estimation of 
the critical temperature T c and the correlation length critical exponent {u)) [33]. To this, 
we Define the following quantities: 



Vi = 4[M 3 ] -3[M 4 ], 


(40) 


V 2 = 2[M 2 ] - [M 4 ], 


(41) 


V 3 ee 3[M 2 ] -2[M 3 ], 


(42) 


V 4 ee (4[M] - [M 4 ])/3, 


(43) 


V 5 ee (3[M] - [M 3 ]/2, 


(44) 


\/ 6 = 2[M]-[M 2 ], 


(45) 



where M = Nm is the total magnetization of the system and 

[Ml ee In (46) 

From Eq. (37) it is easy to show that 

Vj « (l/i/) InL + V^tL 1 /"), (47) 

for j = 1, 2, • • •, 6. At the critical temperature (t — 0), Vj should be constants, independent 
of the system size L. Using Eq. (47) one can find the slope of quantities V\ to V e (Eq. 40- 
45) versus ln(L) for the region near the critical point. Scanning over the critical region and 
looking for a quantity-independent slope gives us both the critical temperature T c and the 
correlation length exponent v with high precision. Figures (12) and (13) give the examples 
of such an effort for the set of the coupling D/\J\ = 0.2. From these figures, we estimate 
that v = 0.842(2) and T c = 1.198(1). The linear fits to the data in Fig.(12) have been 
obtained by the linear least squares method. 

Once v and T c are determined accurately, we can extract other static critical exponents 
related to the order parameter (/3) and susceptibility (7). The ratio j3/v can be estimated 
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by using the size dependence of the order parameter at the critical point given by Eq.(37). 
Fig. (14) shows the log-log plots of the size dependence of the order parameter corresponding 
to D/\J\ = 0.5 and D/\ J\ = 0.2. From this figure the ratio (3/v can be estimated as the slope 
of the straight lines fitted to the data according to Eq.(37). We then have (3/v = 0.198(8) 
for D/\J\ = 0.5 and (3/v = 0.285(8) for D/\J\ = 0.2. 

Accordingly, from Eq.(38) it is clear that the peak values of the finite-lattice susceptibility 
(x = ((M 2 ) — (M) 2 )/(NT)) and the magnitude of the true susceptibility at T c (the same 
as x w hh (m) = 0) are asymptotically proportional to LPl". Then the slope of straight line 
fitted linearly to the log-log plot of these two quantities versus linear size of the lattices, 
can be calculated to estimate the ratio j/v. In Fig. (15) the finite lattice susceptibility have 
been depicted for D/\J\ = 0.2,0.5, respectively. The slopes of linear lines fitted to these 
data give 7/1/ = 1.39(2) for D/\J\ = 0.5 and 7/1/ = 1.42(2) for D/\J\ = 0.2, where the error 
includes the uncertainty in the slope resulting from uncertainty in our estimate for T c . 

The above procedure has been applied for other values of D/\J\ = 1.0,0.5,0.2 and the 
obtained critical exponents are listed in Table. (II). In this table, the critical exponent a, has 
been calculated using the hyper-scaling relation: 

a = 2-dv, (48) 

in which d = 2. On the other hand the Rushbrook scaling law (a + 2(3 + 7 = 2) is satisfied 
for all set of exponents within the computational errors. For comparison, we have listed 
the corresponding critical exponents of Onsager's solution for 2D-Ising, and also Zamolod- 
chikov's conjecture for the Ising-tricritical point in two-dimensions, which corresponds to a 
2D-0 6 field theory [34]. Zamolodchikov's conjecture is based on conformal field theory and 
has been verified by Monte Carlo simulation [35]. 

One can see from Table. (II) that the critical exponents for D/\J\ = 1.0 are pretty close 
to the 2D-Ising values, then anisotropy magnitude of D/\ J\ = 1.0 is large enough to suppress 
the transverse fluctuations normal to easy-axes directions. Upon decreasing the anisotropy, 
the transverse fluctuations become important and the exponents deviate from Ising values. 
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However, although the exponents u, 7 and a monotonously tend to the the 2D-triciritcal 
values, but the exponent (3 gets farther from it. This discrepancy, might the sign of a new 
universality class, other than 2D-0 6 model. 

At the end, we deal with the dependence of the transition temperature to the anisotropy 
intensity. We have already mentioned the method of obtaining the critical temperature for 
the continuous transitions (D > 0.17). For strongly enough first order transitions whose 
energy histograms are double peaked (D ( 0.15), the finite size transition temperatures 
(T C (L)) , are determined as the temperature at which the two peaks have equal heights. 
Once T C (L) for all lattice sizes is obtained, the transition temperature in thermodynamic 
limit can be extrapolated by the following scaling relation: 



where B is a constant and d = 2. The resulting transition temperatures are listed in 
Table. (I). In Fig. (16), we have plotted the transition temperature versus D in logarithmic 
scale. This linear log-log plot shows a power law relation between these to quantities as: 



This result is in clear contrast with mean-field prediction of a linear dependence of transition 
temperature on the anisotropy intensity D. This scaling behaviour can be explained by a 
simple dimensional analysis. Assuming that both exchange interaction, J, and anisotropy, 
D, are equally important in occurrence phase transition in XY Kagome antiferro magnet. So 
the thermal energy which balances the entropy and internal energy at the transition point, 
must be proportional to a combination of J and D. Accordingly, dimensional analysis 
requires K B T C ~ (\J\D)2, which leads us to T C /\J\ ~ (D/\ J|)5. 




(49) 



T c OC £)°- 501 ( 2 ). 



(50) 



IV. CONCLUSION 



In summary, using mean-field theory and the optimized Monte Carlo simulation based on 
multi-histogram, we investigated the phase transitions of the antiferromagnetic classical XY 
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model on a two dimensional Kagome lattice with the easy-axes single ion anisotropy. In the 
absence of anisotropy this system is highly frustrated and no phase transition is expected to 
occur at finite temperatures, except the Kosterlitz-Thouless transition mentioned in Ref. [24]. 
Turning on the anisotropy, removes the degeneracies of the ground state and so establishes 
a long range order with all-in all-out spin configuration at low temperatures. By increasing 
the temperature, the system exhibits a phase transition from all-in all-out ordered state 
to disordered (paramagnetic) state. According to Monte Carlo results this transition is 
first order for small values of anisotropy, while turns to second order at a tricritical point, 
corresponding to an anisotropy strength in the interval 0.15 < Mr < 0.17. 

Employing finite size scaling theory, we derived the critical exponents for continuous tran- 
sitions and found that the transition is in Ising universality for large values of anisotropy. 
This is because in large D/\J\ limit, the fluctuations perpendicular to easy- axes directions 
are frozen, and so the effective degrees of freedom are spin flips along easy-axes directions, 
such that the order parameter possess the discrete Z 2 symmetry. Decreasing the anisotropy 
magnitude, activates the spin fluctuations perpendicular to the easy-axes directions. In 
principle, the coupling of transverse modes (independent of anisotropy) and also of other 
underlying modes, shown in Fig. (2) and (3), with the all-in all-out state at q — 0, is the 
reason for the deviation of the universality class of transitions from Ising, and is also re- 
sponsible for changing the type of transition to dis-continuous for small values of anisotropy. 
However, obtained critical exponents near the tricritical point, do not coincide with those 
of two-dimensional Ising-tricritical point derived from 2D-0 6 field theory. This suggests the 
possibility of the existence of a new tricritical universality class in two- dimensions. It is not 
surprising, because the critical behaviours in frustrated systems are usually different form 
standard universality classes [36]. In this case, finding such a universality class requires more 
theoretical and numerical investigations. 

We hope that this work will motivate further experimental, computational and analytical 
efforts for deeper understanding of the nature of transitions in geometrically frustrated 
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systems. 
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TABLES 



n / 7 
V/J 




U 


1.0 


0.449(1) 


0.66662(7) 


n r 

0.5 


0.316(1) 


0.66660(9) 


0.2 


0.199(1) 


0.66659(8) 


0.18 


0.189(5) 


0.66653(9) 


0.17 


0.184(b) 


0.66649(9) 


0.15 


0.174(8) 


6664(1"! 


0.14 


0.167(7) 


0.6662(1) 


0.13 


0.162(7) 


0.6661(1) 


0.12 


0.156(8) 


0.6659(1) 


0.1 


0.142(8) 


0.6658(1) 



TABLE I. The critical temperatures and value of U* for 
^ = 1.0, 0.5, 0.2, 0.18, 0.17, 0.15, 0.14, 0.13, 0.12, 0.1. (see the text) 



D/\J\ 


V 


/? 


7 


a 


a + 2/3 + 7 


1 


1.019(2) 


0.15(1) 


1.64(8) 


-0.038(4) 


1.9(1) 


0.5 


0.959(2) 


0.19(1) 


1.52(6) 


0.082(4) 


2.0(1) 


0.2 


0.842(2) 


0.24(2) 


1.18(6) 


0.316(4) 


2.0(1) 


2D-Ising 


1 


1/8 


7/4 


0(log) 


2 


2D-c/> 6 


5/9 


1/24 


37/36 


8/9 


2 



TABLE II. The static critical exponents v, (3, 7 and a for -j = 1.0,0.5,0.2, derived from finite- 
size scaling. In the last column the Rushbrook's scaling law is computed. The last two rows are 
listed the corresponding exact critical exponent of 2D-Ising model and two-dimensional Ising-tri- 
criticl point, respectively. 
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FIGURES 
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FIG. 2. spectrum of coupling matrix J for D = 0.2 along [10] direction, degeneracies have been 
removed by addition anisotropic term. 
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FIG. 3. spectrum of coupling matrix J for D = 1.0 along [10] direction, degeneracies have been 
removed by addition anisotropic term. 
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FIG. 4. All in-all out configuration in kagome' lattice. 
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FIG. 5. Temperature dependence of Energy per spin for D = 0.2, 0.1. 
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FIG. 6. Temperature dependence of order parameter (magnetization) for D = 0.2,0.1. 
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FIG. 7. Temperature dependence of specific heat for D = 0.2,0.1. 



29 



f I 

300 - 



y 200 - 



100 - 



0- 




FIG. 8. Temperature dependence of susceptibility for D = 0.1, 0.2. 
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FIG. 9. Size dependences of binder's fourth energy cumulant for D = 1.0,0.2,0.15,0.13,0.1. 
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FIG. 10. Three energy histograms for D = 0.2 and size N = 3 x 40 x 40 near the transition 
temperature. 
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FIG. 11. Energy histogram for D = 0.1 and size ./V = 3x40x40 near the transition temperature. 
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FIG. 12. Dependence of quantity Vj (see the text) versus logarithm of L for D = 0.2 at 
T = 0.1989(5). The solid lines represent linear fits to Eq.(47). All straight lines have the same 
slope v = 0.842(2). 
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FIG. 13. Scanning results for the dependence of quantity Vj versus j for D = 0.2. The hori- 
zontal line is drawn at 1/u = 1.187. 
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FIG. 14. Log-Log plot of order parameter for D = 0.5,0.2. The slopes of fitted line gives 
§ = 0.285(7) for D = 0.2 and £ = 0.198(8) and for D = 0.5. 
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FIG. 15. Log-Log plot of finite lattice susceptibility for D = 0.2, 0.5. The slopes of fitted lines 
give I = 1.58(5) for D = 0.2 and £ = 1.40(5) for one D = 0.5. 
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